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Summary 

The purpose of this paper is to describe an accurate, yet 
economical, method for computing temperature effects in laminar 
lubricating films in two dimensions. The procedure presented here is 
a sequel to one presented in Leeds in 1986 that was carried out for 
the one-dimensional case. 

Because of the marked dependence of lubricant viscosity on 
temperature, the effect of viscosity variation both across and 
along a lubricating film can dwarf other deviations from "ideal" 
constant-property lubrication. 


In practice, a thermohydrodynamics program will involve 
simultaneous solution of the film lubrication problem, together 
with heat conduction in a solid, complex structure. The extent of 
computation required makes economy in numerical processing of 
utmost importance. In pursuit of such economy, we here use techni- 
ques similar to those for Gaussian quadrature. We show that, for 
many purposes, the use of just two properly positioned temperatures 
(Lobatto points) characterizes well the transverse temperature 
distribution . 
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X. INTRODUCTION: 


The purpose of this paper is to describe an accurate, yet 
economical, method for computing temperature effects in laminar 
lubricating films. The procedure presented here is a sequel to one 
presented in Leeds in 1986. 1 

Because of the marked dependence of lubricant viscosity on 
temperature, the effect of viscosity variation both across and 
along a lubricating film can dwarf other deviations from "ideal" 
constant-property lubrication. In two recent papers 2>3 Khonsari has 
summarized the growing literature concerned with the present 
subject. Consequently, we shall not undertake a survey here, but 
refer only to those articles used for support or comparisons. 

In practice, a thermohydrodynamics program will involve 
simultaneous solution of the film lubrication problem, together 
with heat conduction in a solid, complex structure. The extent of 
computation required makes economy in numerical processing of 
utmost importance. In pursuit of such economy, we here use techni- 
ques similar to those for Gaussian quadrature. We show that, for 
many purposes, the use of just two properly positioned temperatures 
characterizes well the transverse temperature distribution. 


2. NOMENCLATURE 

Single-underline _ denotes a Legendre coefficient. 
Double-underline _ denotes a vector. 

Certain matrices are defined in place in section 7 . 


A flow vector, defined by eq. 4.03 

B flow vector, defined by eq . 4.04 

C P specific heat, J/kg-C 

h film thickness, m 

k thermal conductivity, J/Csm 

l subscript for lower surface 

m mass flux vector, kg/sm 2 

p pressure, N/m 2 

Pk Legendre polynomial of order k 
t time, sec 

T temperature, deg C 

u velocity of fluid in x-direction, m/s 
ii subscript for ]ipper surface 

v velocity of fluid in y-direction, m/s 
v total velocity vector, e * u + e y v + £2 w 
V velocity vector, e.x u + e>v, m7s 
w velocity of fluid in z-cTirection , m/s 
Wk weight for ordinate at <Tu , Lobatto quadrature 
x longitudinal position along film, m 
y lateral position along film, m 

z position normal to film, measured from midsurface 


2 



$ dissipation, J/sm 3 , defined by eq. 4.03 

p viscosity, Ns/m z 

C fluidity = 1/p, m 2 /Ns 

2z/h, fractional transverse position 
p density, kg/m 3 

3. LOBATTO INTERPOLATION AND QUADRATURE: 

In general terms, the problem to be treated is one of three- 
dimensional heat convection. A numerical solution is effected by 
sampling the velocities, pressure and temperature over a chosen 
grid of points, and interlinking their values by algorithms that 
incorporate appropriate physical laws. Generally speaking, the 
fewer sample points required, the less the computational effort. 
Consequently, we ask: "What are the best cross-film sampling 

positions?" The answer is provided by the theory of orthogonal 
polynomials . 4 • 5 

Figure 1 shows a section of lubricating film, with the normal 
displacement, £ , scaled to be -1 on the lower wall, and +1 on the 
upper. At certain locations, £k , we intend to obtain sample values 
of the flow variables and to deduce therefrom certain intermediate 
cross-film values, derivatives and integrals. For illustrative 
purposes, consider the temperatures Tk = T(£k). These are assumed 
to be known at i = l and <r = -l , and at N intermediate points. 

If the intermediate sample points are equispaced, and a 
polynomial for TU) is passed through them, the excursions of such 
a polynomial between points can become unacceptably large. See, for 
example, Fig. 2, where a tenth-degree polynomial is used to 

approximate the stepfunction T(<f<0) = 1; T(£>0) = 0 by collocation 
at eleven equispaced points. On the other hand, if a tenth-degree 
polynomial is collocated at the endpoints (<T = —1; ^ = 1) and at 
the zeroes (£k) of the Jacobi polynomial ° ) Pa 1 • 1 ( <T ) ("Lobatto 
points") , conformity to the step function is much better. Indeed, 
if the order of the polynomial is increased, equispace interpola- 
tion may fail to converge at all, whereas the latter form of 
interpolation becomes progressively better. 

Not only do the Lobatto points serve well for interpolation 
by high-order polynomials, but they also serve for accurate 
numerical integration 5 . It can be shown that N such internally- 
selected points permit exact numerical integration of a polynomial 
of order 2N+1 over the range -1<£<1. Thus: 


< a > For N interior points, the Lobatto locations (£i. ) are at the 
zeroes of dP N 4 t U)/d<r = N (N+l ) P N - 1 1 • 1 U) , where Pn (<T) is the Legendre 
polynomial of order N. 
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1 

[2.01] JT ( £ ) d£ = ZwkTk 
-1 

The Lobatto locations, £k , and weight factors, Wk , are to be found 
in Abramowitz and Stegujfa . For N=2, we have: 

Location 
-1 

-1/X5 
1/X5 
1 

Note that, in contrast to Gaussian quadrature, Lobatto 's technique 
incorporates endpoint values which, in our case, are known, and 
might as well be used. 

The approach taken in this paper is to form partial differen- 
tial equations in x,y for the Lobatto-point temperatures, with the 
transverse temperature distribution given by a collocated polynomi- 
al. The "fluidity", £ = 1 / p is also collocated to its Lobatto-point 
values , C k = £ (Tk ) . 

Let us turn now to the analytical development which incor- 
porates these ideas. 

3. BASIC THERMOHYDRODYNAMIC EQUATIONS: 

The equations to be used are for laminar lubrication with a 
fluid possessing constant density and constant thermal properties. 
The momentum equation is: 

[4.01] (a/az) (pav/az) = ?P 

the pressure being independent of "z". 

And the energy equation is: 

[4.02] pC P (DT/Dt) = k(a 2 T/az 2 ) + $ with: 

[4.03] $ = pOv/dz) 2 

Here the viscosity, p, may depend markedly on the local temperature 
T. In addition, the conservation of mass (or volume) must be 
satisfied; i . e . , 

[4.04] n*v = 0 

Several investigators 7 ■ 8 have justified the use of these 
equations by showing that the effects of variations in density, 
specific heat, etc., of the lubricant are quite secondary to those 
resulting from changes in its viscosity, which usually varies 
substantially with temperature. 
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The boundary conditions to be satisfied by eqs . [4.01 - 4.03] 
are as follows. The liquid velocities must conform to the upper 
and lower surface motions, and the liquid temperature to the 
surface temperatures. At the film peripheries, incoming liquid 
temperatures correspond to ambient conditions, whereas outgoing 
temperatures derive from the interior where they originate. 

5. VELOCITIES. REYNOLDS EQUATION: 

To compute the velocities, it is, for several reasons, conven- 
ient to use the fluidity expressed as a series in Legendre 
polynomials. Thus: 


[5.01] £ = ZikPu (f > 


As previously stated, the coefficients fji are here determined by 
collocation of the fluidity at the Lobatto points, where its values 
are determined by the Tk . A double integration of eq . [5.01] then 
yields for V = e*u + v . 

<r <r 


[5.02] 
where : 

[5.03] 
and : 


v = v t + a jc at + b m d<r 

" -l -l 

l l 

A = [Vu - Vi - B JV£d<r] / [JCd<r] 

-1 -1 


[5.04] B = (h/2 ) 2 9P 

To obtain the lineal mass flux, eq. [5.02] must be integrated 
again. The result is: 


[5.05] m/p = (Vu + Vl ) (h/2 ) - (h/3) Li A - (2/3) B (Lo + 2£_2 /5 ) 

This result is independent of the number of terms in the series 
[5.01]. When it is inserted into the mass continuity equation: 


[5.06] ah/at + y*(m/p) = 0 

the following generalized Reynolds equation results: 

[5.07] <7 * £ p h 3 <7p = 6(Vu + Vl ) * *7h + 12(ah/at) - 29* ( £ i /£ o ) h ( Vu - Vi ) 
where : 

[5.08] £p = L° +0.4£_2 _ (La) 2 /<3Lo) 

In form, eq. [5.07] differs from the standard Reynolds equation 9 
only through its last term. 
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Now we wish to express our differential equations using the 
fractional gap position as transverse coordinate across the film. 
Thus : 


[5.09] ^ = 2z/h 

where "z" is measured from the midsurface. Taking "x" as a typical 
lateral coordinate and "u" as a typical variable, we note that: 

[5.10] (8u/dx)y,z = {3u/3y)y,£ - < { 3log { h ) /dx) y ( 3u/d<T ) x , y 

The transverse velocity, w, can be found from mass continuity 

via : 

[5.11] aw/az = -(au/ax + av/ay) = -(9*^)* 

But : 


[5.12] (?‘V) 2 = (^-V)^ - <T91og (h) * (aV/ a<T ) 

Substituting [5.12] into [5.11] and integrating, we get: 


z 

[5.13] w = w L - X (<7<r*v - ^yiog(h) ♦ (av/a^) idz 


-h/2 

Integration by parts gives: 


<r 

[5.14] w = W L + V| *7(h/2) + <rV*<7(h/2) - • [ (h/2) XVd<T] 

” -1 


6. TEMPERATURE EQUATION: 

To obtain a convenient differential equation for the tempera- 
ture, we rewrite eq. [4.02] as: 

[ 6 . 01 .] 

aT/at + u(aT/ax)y.^ + v(aT/ay)*, f + (2/h) (aT/a«r)x,y [w -fV*V(h/2) i = 

(4ic/h 2 ) (a 2 T/d£ 2 )x,y + $/(pCp ) 

To treat the cross-velocity term, we note that kinematics gives: 

[6.02] W L + V* ^ (h/2 ) = -a(h/2)/at 

Substituting [5.14] and [6.02] into [6.01], we get: 

<: 

[6.03] aT/at + v^r - (l/h) (bt/ 3£) [(i+<r) (ah/at) + 9*h/vd<r] = 

" -l 

(4K/h 2 ) (a 2 T/a<r 2 ) + $/(pC P ) 
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All partial derivatives in the foregoing equation are in x,y,<T,t 
space. There is one such partial differential equation (6.03] in 
x,y,t for each Tk . 

7. NUMERICAL PROCEDURES: 

It is convenient to treat the Tk (xi , yj ) , £k (xi , yj ) as vector 
matrices, and to use matrix notation. We have: 


[7. 

.01] 

T(<r) 

= ZTn 

Pu U) 

[7. 

.02] 

£ (<r) 

3 

II 

Pu (C) 


Therefore : 

[7.02] Tk = ZTn PuUk) or: Tk = ZCk ■> Tk or: T = CT and: T = C' l T 
Also, then: 

[7.03] £ = C£_ and: L = C- ' £ 

Successive differentiation of [7.01] enables us to write: 

[7.03] aT/a^ = DT and: a z T/a<T z = ET 

Successive integration of [7.02] permits us to write 

[7.04] V = Vl + AF£ + BG£ and: 

sf 

[7.05] JVd<r = Vi(l + <r) + AR£ + BS£ 

-1 = 

To preserve numerical stability, backward differences are used 
for the convective terms . Thus , at the point (i,j) we write: 

[7.06] u(aT/3x) = abs(ui.j) [Ti.j - Tt-gx.j]; gx = sgn(ui,j) 

with a similar expression for v(aT/ay). 

Now form the diagonal matrix, Ea , from the diagonal terms of 
E, and let E 0 = E - E,t . With these definitions, the matrix finite- 
difference equation corresponding to [6.03] is: 

[7.07] [1/At + abs(u)/Ax + abs(v)/Ay - ( 4 k / h z ) Ea ] T new i , j = 

Ti.j/At + { abs (u) Ti -gx , j 1 /Ax + 1 abs ( v ) Ti , j - gy I / Ay + (4K/h 2 )E 0 Ti.j + 

( 1/h) (yhjVd^)D Ti , j + $/pC P 
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Here T new denotes T(t+At), u, v, $ and T belong to the same trans- 
verse row (£k ) , and the integral of velocity is from the lower wall 
to <Tk - 


The explicit form of the difference equation [7.07] is stable 
for a large At, such as was used to generate the steady-state 
results presented in this report. 

8. COMPUTATIONAL RESULTS: 

Validation of the computer program was attempted in various 
ways. Comparisons were made with the one-dimensional calculations 
of Hunter and Zienkiewicz 1 0 , Dowson and Hudson 7 and of Hahn and 
Kettleborough 8 , who all used the physical properties tabulated 
below . 

p = 0.13885 exp 1 -y (T-Tanb i en t ) I , Ns/m 2 ; P = 0.045 
pC P = 1.7577 J / m :) C 
k = 7 . 306 E-08 m 2 /s 

Ta nbtent = 0 (used as reference, only) 

These same properties are used for all examples of this report, 
except that p * 0, instead of 0.045, for the constant-property 
calculations . 

Dowson and Hudson chose the following bearing characteristics. 
L = 0.18288 m; W infinite 
Uu = 31.96 m/s; Ul = 0 

hi = 1.8288 E-04 m; h 2 = 0.9144 E-04 m 

Figures 3 and 4 compare our results with theirs (b) . Agreement is 
certainly satisfactory, though not perfect. Dowson and Hudson used 
Ax = L/20 and a constant Ay =h 2 /20, whereas we employed Ax = L/30 
and N = 8. In these comparisons, as in ot he rs made, the reasons for 
the observed small discrepancies are difficult to assign. There- 
fore, it was decided to write the present program for arbitrary N, 
and to use N = 8 as a standard against which to test more approxi- 
mate, but faster versions employing N = 2 or 3 . To assist others 


(b > All dotted curves, except those in Fig. 2, were collocated at 
the Lobatto points, and interpolated with an auxiliary plotting 
program; i.e. the dotted curves are not everywhere in complete 
agreement with the corresponding Legendre polynomial expansions. 
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who may wish to make comparisons with their own programs, we list 
some results for the above problem in the Appendix. 

In Fig. 3, the great reduction in load capacity due to 
lubricant warmup is manifest. Of course, in engineering design 
practice some allowance for this effect is made by assigning some 
lowered constant viscosity corresponding to an estimated tempera- 
ture rise. 

Figure 4 shows the convergence with N of results obtained 
using the present program. There is no perceptible difference 
between the results for N = 5 and those for N = 8. Somewhat 
fortuitously, the results for N = 2 also almost coincide with the 
"true" curve. Figure 5 shows the convergence of the temperature 
distributions. The Lobatto-point results for N = 2 and N = 3 are 
compared with the curves for N = 8 at the bearing exit, and at the 
halfway point. 

Reverse convection has been a source of difficulty for some 
investigators. Accordingly, we show in Figs. 6 and 7 some calcula- 
tions for a high film-thickness ratio of 4. Note the rapid 

variation of temperature in the bearing inlet incoming 

temperatures are taken at an entrance value of 0, whereas tempera- 
tures in the backflow region originate from a region where viscous 
heating has occurred. As shown in Fig. 1, Lobatto interpolation is 
particularly suited to cope with such variation. The two-tempera- 
ture version (N - 2) of our program performs surprisingly well. 

All of the foregoing calculations are for a one-dimensional 
bearing, run with the program, however, as a two-dimensional wide 
bearing. The same cases were also solved earlier by a one-dimen- 
sional procedure 1 embodying Lobatto-point methods. But when 
extended to two-dimensional problems, that earlier procedure proved 
to be persistently unstable. To deal successfully with the added 
dimension, the present method was devised. Figures 8, 9 and 10 show 
some sample "true" two-dimensional results obtained with a square 
slider bearing. 

We intend to continue this investigation by engaging in some 
parametric studies and by coupling the new technique to: 

inlet-groove temperature distribution 
cavitation 

heat conduction in the bearing body 
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APPENDIX: 


PRESSURE DISTRIBUTION FOR DOWSON-HUDSON PROBLEM 


£ ZL 

oressure, N/: 

0 

0 

. 1 

3545479 

.2 

6434485 

. 3 

8755968 

. 4 

1 . 051984E+07 

. 5 

1 . 167347E+07 

.6 

1 . 2097 84E+07 

.7 

1 . 158809E+07 

. 8 

9815331 

.9 

6259130 

1 

0 


TEMPERATURES OBTAINED FOR DOWSON-HUDSON PROBLEM 


£ 

-1 

-.9195339 

-.7387739 

-.4779249 

-.165279 

.165279 

.4779249 

.7387739 

.9195339 

1 


x/L 


0 

0 . 2 

0 . 4 

0.6 

0 . 8 

1 . 0 

0 

0 

0 

0 

0 

0 

0 

1 . 883201 

3.275796 

4 .859227 

6.956678 

9 .435757 

0 

5 . 325756 

S . 360686 

13 . 81515 

19.56049 

26 .13743 

0 

7.951722 

14 . 41125 

21.01087 

28.58318 

36 .01018 

0 

8 . 313656 

15 . 44364 

21 . 92048 

27.68786 

31 . 61222 

0 

7 .881326 

14 . 21572 

19 . 14167 

22-08036 

22.90305 

0 

7 . 5885 

12 . 91233 

15.74303 

16.07165 

15.18215 

0 

6.967654 

10 . 07563 

10.47426 

9.308805 

8 . 233334 

0 

3 . 865244 

4 . 300633 

3 . 853627 

3.094525 

2.774959 

0 

0 

0 

0 

0 

0 
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Figure 4. Temperature Contours (T/311) for 
Infinite-Width Slider Bearing. 
Dowson- Hudson ( ) vs . N=8 (X) 

Variable Viscosity, hi /ha = 2 
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Figure 5. Convergence of Pressure Profiles with N 
Variable Viscosity, Infinite-Width Slider, hi /h 2 = 2 

N=8 (...); N= 5 (0); N=3 (*); N=2 (X) 
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TEMPERATURE 



- 6,0 - 2.0 2.0 6.0 10 . 0 * E-i 

POSITION IN GAP 


Figure 6. Temperature Profiles at Exit and Halfway 
Infinite-Width Slider, Variable Viscosity, hi /h 2 = 2 

N=8 N=3 (*) ; N=2 (X) 




2.0 <.0 6.0 6.0 10 . 0 * 8—1 

LONGITUDINAL POSITION 


Figure 7. Pressure Profile for Infinite-Width Slider 
Variable Viscosity, hi /hi = 4 
N=8 (...); N=2 (X) 
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Figure 8. Velocity Profiles at Exit, Halfway & Entrance 
for Infinite-Width Slider, hi /hi = 4 
N=8 (...); N=2 (X) 
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Figure 11. Center Temperature Protiles for Square Slider 

Variable Viscosity, hi /hz 
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